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Abstract 

Let A(x) = Aq + x\A\ + • • • + x n A n be a linear matrix, or pencil, generated by 
given symmetric matrices Aq,A\,... ,A n of size m with rational entries. The set 
of real vectors x such that the pencil is positive semidefinite is a convex semi- 
algebraic set called spectrahedron, described by a linear matrix inequality (LMI). 

We design an exact algorithm that, up to genericity assumptions on the input 
matrices, computes an exact algebraic representation of at least one point in the 
spectrahedron, or decides that it is empty. The algorithm does not assume the 
existence of an interior point, and the computed point minimizes the rank of the 
pencil on the spectrahedron. The degree d of the algebraic representation of the 
point coincides experimentally with the algebraic degree of a generic semidefinite 
program associated to the pencil. We provide explicit bounds for the complexity 
of our algorithm, proving that the maximum number of arithmetic operations that 
are performed is essentially quadratic in a multilinear Bezout bound of d. When 
m (resp. n) is fixed, such a bound, and hence the complexity, is polynomial in 
n (resp. m). We conclude by providing results of experiments showing practical 
improvements with respect to state-of-the-art computer algebra algorithms. 

Keywords: linear matrix inequalities, semidefinite programming, computer algebra al¬ 
gorithms, symbolic computation, polynomial optimization. 


1 Introduction 


Let S m (Q) be the vector space of m x m symmetric matrices with entries in Q. Let 
A 0 , Ai, ..., A n G § m (Q) and A(x) = Aq+x x A x +- ■ ■+x n A n be the associated linear matrix , 
x = (xi,..., x n ) G M n . We also denote the tuple (Aq, A x , ..., A n ) by A G §” +1 (Q). For 
x G M n , A(x) is symmetric, with real entries, and hence its eigenvalues are real numbers. 
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The central object of this paper is the subset of i 6 1" such that the eigenvalues of A(x) 
are all nonnegative, that is the spectrahedron 

y = {i e R n : A{x) y o}. 

Here y 0 means “positive semidefinite” and A(x) y 0 is called a linear matrix inequality 
(LMI). The set 5? is convex closed basic semi-algebraic. This paper addresses the following 
decision problem for the spectrahedron JAL 

Problem 1 (Feasibility of semidefinite programming) Compute an exact algebraic 
representation of at least one point in SA, or decide that 5? is empty. 

We present a probabilistic exact algorithm for solving Problem (JT]). The algorithm de¬ 
pends on some assumptions on input data that are specified later. If 5? is not empty, 
the expected output is a rational parametrization (see e.g. [55]) of a finite set Z C C n 
meeting 5? in at least one point x*. This is given by a vector (q 0 , ..., q n + i) C ’L\t] n+2 and 
a linear form A = Aiaq + • —I- X n x n e Q[t] such that deg(g n+ i) = jj Z, deg [qf) < deg(g n+ i) 
for 0 < i < n, gcd(g n+ i, q 0 ) = 1 and Z coincides with the set 

{(xi,... ,x n ) E C n | t = AiXi H-h X n x n ,q n+1 (t ) = 0 ,Xi = qi{t)/q Q {t)} (1) 

A few remarks on this representation are in order. Usually, g 0 is taken as for a 

better control on the size of the coefficients m Theorem 1] - see also the introductory 
discussion of that theorem in H3j. More precisely, if D bounds the degrees of a finite family 
of polynomials in Z[xi,..., x n ] defining Z and h bounds the bit size of their coefficients 
and those of A, then the coefficients of a rational parametrization encoding Z have bit 
size bounded by hD n . 

This is to be compared with polynomial parametrizations where the rational fractions 
qi/qo are replaced by polynomials pp they are obtained by inverting q 0 w.r.t. q n+1 using 
the extended Euclidean algorithm. That leads to polynomials pi with bit size bounded 
by hD 2n . 

In order to compute such representations, the usual and efficient strategy is to compute 
first the image of such representations in a prime field and next use a Newton-Hensel 
lifting to recover the integers. According to [231 Lemma 4] and the above bounds, the 
cost of lifting integers is log-linear in the output size. Since in the case of polynomial 
parametrizations, the output size may be D n times larger than in the case of rational 
parametrizations, rational parametrizations are easier to compute. In addition, observe 
from [1451 Lemma 3.4 and Theorem 3.12] that isolating boxes for the real points in Z from 
rational or polynomial representations have the same bit complexity (e.g. cubic in the 
degree of q and log-linear in the maximum bit size of the coefficients in the parametriza¬ 
tion). 

As an outcome of designing our algorithm, we also compute the minimum rank attained 
by the pencil on the spectrahedron. Moreover, since the points in Z are in one-to- 
one correspondence with the roots of q n +i, from this representation the coordinates of 
the feasible point x* G 5? can be computed with arbitrary precision by isolating the 
corresponding solution t* of the univariate equation q n+ \(t) =0. If 5? is empty, the 
expected output is the empty list. 
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1.1 Motivations 


Semidefinite programming models a large number of problems in the applications [53, 
m i- This includes one of the most important questions in computational algebraic 
geometry, that is the general polynomial optimization problem. Indeed, Lasserre [12 ] 
proved that the problem of minimizing a polynomial function over a semi-algebraic set can 
be relaxed to a sequence of primal-dual semidefinite programs called LMI relaxations, and 
that under mild assumptions the sequence of solutions converge to the original minimum. 
Generically, solving a non-convex polynomial optimization problem amounts to solving a 
finite-dimensional convex semidefinite programming problem mi- Numerical algorithms 
following this approach are available and, typically, guarantees of their convergence are 
related to the feasibility (or strict feasibility) of the LMI relaxations. It is, in general, a 
challenge to obtain exact algorithms for deciding whether the feasible set of a semidefinite 
programming (SDP) problem 


min > CiXi s.t. A(x) >z 0 (2) 

iei" ^ 

2=1 

is empty or not. The feasible set of the SDP ([2]) is defined by an LMI and hence it is a 
spectrahedron. Problem ([I]) amounts to solving the feasibility problem for semidefinite 
programming, in exact arithmetic: given a Q—definable semidefinite program as in (J2J) 
(that is, we suppose that the coefficients of A(x) have rational entries), decide whether the 
feasible set 5? = {x G M n : A(x) >z 0} is empty or not, and compute exactly at least one 
feasible point. We would like to emphasize the fact that we do not assume the existence 
of an interior point in 5?. Quite the opposite, we are especially interested in degenerate 
cases for which the maximal rank achieved by the pencil A(x) in SA is small. 

This work is a first step towards an exact approach to se mi de fini te programming. In 
particular, a natural perspective of this work is to design exact algorithms for deciding 
whether the minimum in ([2]) is attained or not, and for computing such a minimum in the 
affirmative case. While the number of iterations performed by the ellipsoid algorithm ra 
to compute the approximation of a solution of (J2j) is polynomial in the number of vari¬ 
ables, once the accuracy is fixed, no analogous results for exact algorithms are available. 
Moreover, since the intrinsic complexity of the optimization problem ([2]) is related to its 
algebraic degree S as computed in |3H1 23], the paramount goal is to design algorithms 
whose runtime is polynomial in 5. The algorithm of this paper shows experimentally such 
an optimal behavior with respect to 5. 

We finally recall that solving LMIs is a basic subroutine of computer algorithms in systems 
control and optimization, especially in linear systems robust control [21133], but also for 
the analysis or synthesis of nonlinear dynamical systems [67], or in nonlinear optimal 
control with polynomial data [23111]. 

1.2 Contribution and outline 

We design a computer algebra algorithm for solving the feasibility problem of semidefinite 
programming, that is Problem ([T]), in exact arithmetic. Let us clarify that we do not 
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claim that an exact algorithm can be competitive with a numerical algorithm in terms of 
admissible size of input problems: indeed, SDP solvers based on interior-point methods 
mm can nowadays handle inputs with a high number of variables that are out of reach 
for our algorithms. Our contribution can be summarized as follows: 

1. we show that the geometry of spectrahedra understood as semi-algebraic sets with 
determinantal structure can be exploited to design dedicated computer algebra al¬ 
gorithms; 

2. we give explicit complexity and output-degree upper bounds for computer algebra 
algorithms solving exactly the feasibility problem of semidehnite programming; our 
algorithm is probabilistic and works under assumptions on the input, which are 
generically satisfied; 

3. we provide results of practical experiments showing the gain in terms of computa¬ 
tional timings of our contribution with respect to the state of the art in computer 
algebra; 

4. remarkably, our algorithm does not assume that the input spectrahedron y = {x G 
M” : A(x) >z 0} is full-dimensional, and hence it can tackle also examples with empty 
interior. 

The main idea is to exploit the relation between the geometry of spectrahedra, and that 
of the determinantal varieties associated to the input symmetric pencil A(x). Let us 
introduce, for r = 0,..., m — 1, the algebraic sets 

V r = {i 6 C“ : rank A(x) < r}. 

These define a nested sequence T> 0 C T>\ C • • • C D m _ i. The dimension of D r for generic 
linear matrices A is known, and equals n — 1 ) (see Lemma HD- The Euclidean 

boundary dSA of y is included in the real trace of the last algebraic set of the sequence: 
dy C V m _i fi M ra . In particular, for x G dy , the matrix A(x) is singular and one could 
ask which elements of the real nested sequence V 0 flM n C • • • C V m _ i flM n intersect dSA . 

Notation 1 If y = { x E MT : A(x) >z 0} is not empty, we define the integer r(A) = 
min {rank A(x) : x E y} . 

When y is not empty, r(A) equals the minimum integer r such that T> r fl M n intersects 
y. We present our first main result, which states that y contains at least one of the 
connected components of the real algebraic set 'D^a) Hi". We denote by §(^ +1 (Q) = 
§ m (Q) x • • • x S m (Q) the (n + 1)—fold Cartesian product of § m (Q). 

Theorem 2 (Smallest rank on a spectrahedron) Suppose that y 0. Let C be a 

connected component of TV(a) H M n such that C fl y 0. Then C C y and hence 
C C (X\-(t4) \ T^r(A)-l) H M n . 
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We give a proof of Theorem [2] in Section [2] From Theorem [2J we deduce the following 
mutually exclusive conditions on the input symmetric linear pencil A: either 5? = 0, or 
y contains one connected component C of XVpl) HR". Consequently, an exact algorithm 
whose output is one point in the component C C y D X> r (A) would be sufficient for our 
goal. Motivated by this fact, we design in Section [3721 an exact algorithm computing one 
point in each connected component of T> r D M n , for r G {0,..., m — 1}. 

The strategy to compute sample points in V r D R n is to build an algebraic set V r C 
(£n+m(m-r) w } lose projection on the first n variables is contained in T> r . This set is defined 
by the incidence bilinear relation A(x)Y(y) = 0 where Y(y ) is a full-rank m x (m — r) 
linear matrix whose columns generate the kernel of A{x) (cf. Section f3 .11) . Unlike T> ri 
the incidence variety V r , up to genericity conditions on the matrices A 0 , Ai,.. ., A n , turns 
to be generically smooth and equidimensional. The next theorem presents a complexity 
result for an exact algorithm solving Problem dTJ) under these genericity assumptions. 


Theorem 3 (Exact algorithm for LMI) Suppose that for 0 < r < m — 1, the inci¬ 
dence variety V r is smooth and equidimensional and that its defining polynomial system 
generates a radical ideal. Suppose that V r has the expected dimension n — ( m_ J +1 ). There 
is a probabilistic algorithm that takes A as input and returns: 


1. either the empty list, if and only if SA = 0, or 

2. a vector x* such that A(x*) = 0, if and only if the linear system A(x) = 0 has a 
solution, or 

3. a rational parametrization q = (q 0 ,... ,q n +i) C Z[t] such that there exists t* G M 
with q n+ i(t*) = 0 and: 

• A(q 1 (t*)/q 0 (t *),..., q n (t*)/q 0 (t*)) h 0 and 

• rank Aiqfit*)/q 0 (t *),..., q n (t*)/q 0 {t*)) = r(A). 


The number of arithmetic operations performed are in 
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If T/ 0, the degree of q is in O \ +n ) ^ . 


An important aspect of our contribution can be read from the complexity and degree 
bounds in Theorem [3j indeed, remark that when m is fixed, both the output degree 
and the complexity of the algorithm are polynomial functions of n. Viceversa, by the 
constraint n > ( m b[ +1 ) (given by Lemma HJ), one also easily deduces that when n is fixed, 
the complexity is polynomial in m. 

The algorithm of Theorem [3] is described in Sectional Its probabilistic nature comes from 
random changes of variables performed during the procedure, allowing to put the sets T> r 
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in generic position. We prove that for generic choices of parameters the output of the 
algorithm is correct. 

A complexity analysis is performed in Section [5] Bounds in Theorem [3] are explicit ex¬ 
pressions involving m and n. These are computed by exploiting the multilinearity of 
intermediate polynomial systems generated during the procedure, and are not sharp in 
general. By experiments on randomly generated symmetric pencils, reported in Section [6l 
we observe that the output degree coincides with the algebraic degree of generic semidef- 
inite programs, that is with data given in [48J Table 2]: this evidences the optimality of 
our approach. We did not succeed in proving exact formulas for such degrees. 

1.3 Related works 

On input A, 5? can be defined by m polynomial inequalities in Qfaq,... ,x n ] of degree 
< m (see e.g. m)- As far as we know, the state-of-the-art for designing algorithms 
deciding the emptiness of 5? consists only of algorithms for deciding the emptiness of 
general semi-algebraic sets; our contribution being the first attempt to exploit structural 
properties of the problem, e.g. through the smallest rank property (Theorem [2]). 

A first algorithmic solution to deciding the emptiness of general semi-algebraic sets is given 
by Cylindrical Algebraic Decomposition algorithm D3; however its runtime is doubly 
exponential in the number n of variables. The first singly exponential algorithm is given 
in |26j . and has led to a series of works (see e.g. [54, 31]) culminating with algorithms 
designed in [5] based on the so-called critical points method. This method is based on the 
general idea which consists in computing minimizers/maximizers of a well-chosen function 
reaching its extrema in each connected component of the set under study. Applying 
[5] to problem ([T]) requires m 0(jL ^ bit operations. Note that our technique for dealing 
with sets V r is based on the idea underlying the critical point method. Also, in the 
arithmetic complexity model, our complexity estimates are more precise (the complexity 
constant in the exponent is known) and better. This technique is also related to algorithms 
based on polar varieties for grabbing sample points in semi-algebraic sets; see for example 
[21 [31 [581 E9] and its application to polynomial optimization E3- 

To get a purely algebraic certificate of emptiness for 5 ?', one could use the classical ap¬ 
proach by Positivstcllensatz g3l E21 iU- As a snake biting its tail, this would lead to a 
family, or hierarchy, of semidehnite programs [42]. Bounds for the degree of Positivstel- 
lensatz certificates are exponential in the number of variables and have been computed 
in [65J for Schmudgen’s, and in [49] for Putinar’s formulation. In the recent remarkable 
result [44], a uniform 5—fold exponential bound for the degree of the Hilbert 17th prob¬ 
lem is provided. In 0H, an emptiness certificate dedicated to the spectrahedral case, by 
means of special quadratic modules associated to these sets, is obtained. 

All the above algorithms do not exploit the particular structure of spectrahedra under¬ 
stood as determinantal semi-algebraic sets. In [40], the authors showed that deciding 
emptiness of 5? can be done in time 0(nm 4 ) + )) ; that is in polynomial time 

in m (resp. linear time in n) if n (resp. m ) is fixed. The main drawback of this algorithm 
is that it is based on general procedures for quantifier elimination, and hence it does not 
lead to efficient practical implementations. Note also that the complexity constant in the 
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exponent is still unknown. 

Also, in [28], a version of EH dedicated to spectrahedra exploiting some of their structural 
properties, decides whether a linear matrix inequality A(x ) Y 0 has a rational solution, 
that is whether 5? contains a point with coordinates in Q. Remark that such an algorithm 
is not sufficient to solve our problem, since, in some degenerate but interesting cases, 5? 
is not empty but does not contain rational points: in Section 16.21 we will illustrate the 
application of our algorithm to one of these examples. 

As suggested by the smallest rank property, determinantal structures play an important 
role in our algorithm. This structure has been recently exploited in [18] and es for the fast 
computation of Grobner bases of zero-dimensional determinantal ideals and computing 
zero-dimensional critical loci of maps restricted to varieties in the generic case. 

Exploiting determinantal structures for determinantal situations remained challenging for 
a long time. In [31] we designed a dedicated algorithm for computing sample points in 
the real solution set of the determinant of a square linear matrix. This has been extended 
in [36] to real algebraic sets defined by rank constraints on a linear matrix. Observe that 
this problem looks similar to the ones we consider thanks to the smallest rank property. 
As in this paper, the traditional strategy consists in studying incidence varieties for which 
smoothness and regularity properties are proved under some genericity assumptions on 
the input linear matrix. 

Hence, in the case of symmetric matrices, these results cannot be used anymore. Because 
of the structure of the matrix, the system defining the incidence variety involves too 
many equations; some of them being redundant. Hence, these redundancies need to be 
eliminated to characterize critical points on incidence varieties in a convenient way. In 
the case of Hankel matrices, the special structure of their kernel provides an efficient 
way to do that. This case study is done in [35]. Yet, the problem of eliminating these 
redundancies remained unsolved in the general symmetric case and this is what we do in 
Section [3] which is the starting point of the design of our dedicated algorithm. 

1.4 Basic notation 

We refer to 0 H3] [30] H6] for the algebraic-geometric background of this paper. We recall 
below some basic definitions and notation. We denote by M P)(? (Q) the space of p x q 
rational matrices, and GL n (C) the set ofnxn non-singular matrices. The transpose of 
M G M p ^(Q) is M J . The cardinality of a finite set T (resp. the number of entries of a 
vector v) are denoted by j \T (resp. j}u). 

Let x — (aq,..., x n ). A vector / = (/ l5 ..., f s ) C Q[x] is a polynomial system, (/) C Q[x] 
its ideal and Z((f)) = {i£ C n : fi(x) — 0, i — 1,..., s} the associated algebraic set. Sets 
Z((f)) define the collection of closed sets of the Zariski topology of C n . The intersection 
of a Zariski closed and a Zariski open set is called a locally closed set. For M G GL n (C) 
and Z C C n , let M~ X Z = {a: G C n : M x G Z}. With I(S) we denote the ideal of 
polynomials vanshing on S cC" 

Let / = (/i,..., f s ) C Q[x]. Its Jacobian matrix is denoted by Df = (dfi/dxj)^. An 
algebraic set Z C €" is irreducible if Z = Zi U Z 2 where Zi,Z 2 are algebraic sets, 
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implies that either Z — Z\ or Z — Z^. Any algebraic set is the finite union of irreducible 
algebraic sets, called its irreducible components. The codimension c of an irreducible 
algebraic set Z C C n is the maximum rank of Df on Z , where I(Z) = (/). Its dimension 
is n — c. If all the irreducible components of Z have the same dimension, we say that Z is 
equidimensional. The dimension of an algebraic set Z is the maximum of the dimensions of 
its irreducible components, and it is denoted by dim Z. The degree of an equidimensional 
algebraic set Z of codimension c is the maximum cardinality of finite intersections Z fl £ 
where £ is a linear space of dimension c. The degree of an algebraic set is the sum of the 
degrees of its equidimensional components. 

Let Z C C" be equidimensional of codimension c, and let I(Z) = (/r,...,/ s ). The 
singular locus of Z, denoted by sing (Z), is the algebraic set defined by / = (/i ,..., f s ) 
and by all c x c minors of Df. If sing (Z) = 0 we say that Z is smooth, otherwise 
singular. The points in sing (Z) are called singular, while points in reg (Z) = Z \ sing (Z) 
are called regular. Let Z C C n be smooth and equidimensional of codimension c, and 
let I(Z) = (/!,..., f s ). Let g : C n —> C m be an algebraic map. The set of critical 
points of the restriction of g to Z is the algebraic set denote by crit ( g , Z) and defined by 
/ = (/i,..., f s ) and by all c + m minors of the Jacobian matrix D(f,g). The points in 
g(crit ( g,Z )) are called critical values, while points in C m \ g (crit ( g,Z )) are called the 
regular values, of the restriction of g to Z. 

2 The smallest rank on a spectrahedron 

We prove Theorem [2], which relates the geometry of linear matrix inequalities to the 
rank stratification of the defining symmetric pencil. We believe that the statement of this 
theorem is known to the community of researchers working on real algebraic geometry and 
semidefinite optimization; however, we did not find any explicit reference in the literature. 

Proof of Theorem [2} By assumption, the rank of A(x ) on y is greater or equal 
than r(A). We consider the vector function e= (ei,...,e m ) : R n —> R m where 
ei(x) < ... < e m (x) are the ordered eigenvalues of A(x). Let C be a connected component 
of V r (A) n R n such that C fl y ^ 0, and let x G C fl y. One has rank A(x) = r(A) and 
ei(x) = ... = e m _ r p 4 )(a;) = 0 < e m _ r (A)+i{%) < ... < e m (x). Suppose ad absurdum that 
there exists y G C such that y f AA. In particular, one eigenvalue of A(y) is strictly 
negative. 

Let g : [0,1] —> C be a continuous semi-algebraic map such that g(0) = x and g( 1) = y. 
This map exists since C is a connected component of a real algebraic set. The image 
g([0,1]) is compact and semi-algebraic. Let 

T={te [0,1] : g{t) G = ^([O, l])ny). 

Since g is continuous, T C [0,1] is closed. So it is a finite union of closed intervals. 
Since 0 G T (g(0) = x G 5A) there exists t 0 G [0,1] and N G N such that [0,0] £ T 
and for all p > N, t 0 + ^ ^ T. One gets that g(t 0 ) = x G 5? and that for all p > N, 
g[to + = x p f: y. By definition, x,x p G C C TVpi) O R” for all p > N, and since 

x G y, we get rankA(x) = r(A) and rank A(x p ) < r(A) for all p > N. We also get that 


rank A{g(t)) = r(A) for all t £ [0, t 0 ]. We finally have x p —» x when p —> +oo, since g is 
continuous. There exists a map 


<p: |p£N:p>iV|—* Z 

which assigns to p the index of eigenvalue-function among e\,... ,e m corresponding to the 
maximum strictly negative eigenvalue of A{x p ), if it exists, otherwise it assigns 0. Remark 
that since rankH(x p ) < r(A) for all p, then <p(p) < r(A) for all p. In other words, the 
eigenvalues of A(x p ) satisfy 


^ ^ ®<p(p)("R>) ^0 ^<f(p)+l ip^p) ■ ■ ■ W(p)+m—r(,4) i^p) 

0 — ^[p(p)+m—r(A)+l(%p) — • • • — 

for p > N. Since the sequence {(fi(p)} p >N is bounded, up to taking a subsequence, it 
admits at least a limit point by the Bolzano-Weierstrass Theorem jT, Th. 3.4.8], this point 
is an integer, and j t-)- <p(j) is constant for large j. Suppose that there exists a limit point 
£ > 0, and let {pj}jm such that <p(pj) —> £ and that for j > N', j (->■ <p(pj) is constant. 
Thus, 0 = e£ + i(x Pj ) = ... = e£ +m _ r (A)(x Pj ) for all j > N'. Since x Pj —> x, and since 
ei,... ,e m are continuous functions, we obtain that £ = 0 is the unique limit point of <p, 
hence <p converges to 0. Hence <p = 0 for large p. This contradicts the fact that x p ^ 5? 
for large p. 

We conclude that the set C\S? is empty, that is C C 5?. By the minimality of the integer 
r{A) in {rank A(x) : x £ one deduces that C C (V r {A) \ T^r(A)~ i) H R n . □ 


3 Algorithm 

Our algorithm is called SolveLMI, and it is presented in Section [3731 Before, we describe in 
Section [3721 its main subroutine LowRankSym. which is of recursive nature and computes 
one point per connected component of the real algebraic set T> r D M n . We start, in the 
next section, with some preliminaries. 

3.1 Preliminaries 

Expected dimension of low rank loci 

We first recall a known fact about the dimension of V r , when A is a generic symmetric 
pencil. 

Lemma 4 There exists a non-empty Zariski open subset sA C §^ +1 (C) such that, if 
A £ 0§” +1 (Q), for all r — 0,..., m — 1 , the set V r is either empty or it has dimension 


Proof : The proof is classical and can be found e.g. in [H Prop. 3.1]. □ 
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Incidence varieties 


Let A(x) be a symmetric m x m linear matrix, and let 0 < r < m — 1. Let y = 
( Hi,j)i<i<m,i<j<m-r be unknowns. Below, we build an algebraic set whose projection on 
the x— space is contained in T> r . Let 


Y(v) 


^ 2/l,l ' ' ' 2/l, m—r ^ 

y Vm,l ' ' ' 2/m, m—r J 


and let t = {zi,..., i m - r } C {1,..., m}, with j}z = m — r. We denote by Y b the (m — r) x 
(m — r) sub-matrix of Y{y) obtained by isolating the rows indexed by t. There are ( m ) 
such matrices. We define the set 


V r (A, t) = {(x, y) £ C n x C m ^ : A(x)Y(y) = 0, Y t - I m _ r = 0}. 

We denote by f(A,t), or simply by /, when there is no ambiguity on t, the polynomial 
system defining V r (^4, 1 ). We often consider linear changes of variables x: for M e GL n (C), 
f(A o M, i) denotes the entries of A(M x)Y(y) and Y L — I m _ r , and by V r (A o M, t) its zero 
set. We also denote by U b £ M m _ rjm (Q) the full rank matrix whose entries are in {0,1}, 
and such that U,Y(y) = Y t . By simplicity we call U L the boolean matrix with multi-index 
l. 

By definition, the projection of V r {A, l) on the first n variables is contained in T> r . We 
remark the similarity between the relation A(x)Y(y) = 0 and the so-called complemen¬ 
tarity condition for a couple of primal-dual se mi def ini te programs, see for example [JSJ 
Th. 3]. The difference in our model is that the special size of Y(y) and the affine constraint 
Y l = Im-r force a rank condition on Y(y) and hence on A{x). 


Eliminating redundancies 

The system f(A,t) contains redundancies induced by polynomial relations between its 
generators. These relations can be explicitly eliminated in order to obtain a minimal 
polynomial system defining V r . 

Lemma 5 Let M £ GL n (C). Let t C {1,..., m), with ft = m — r. Let A £ Sjy +1 (Q), 
and f £ Q[x,y] m( ' m ~ r ' >+( ' m ~ r ) be the polynomial system defined in Section 1X71 Then we 
can explicitly construct a subsystem f re d C f of length m(m — r) + (”'~ 2 r+1 ) such that 

(fre d ) = (f)- 

Proof : In order to simplify notations and without loss of generality we suppose M = I n 
and l — {1,..., m — r}. We substitute Y t = I m _ r in A(x)Y(y), and we denote by g it j the 
(i, j)— th entry of the resulting matrix. We denote by f re d the following system: 

fred, i j , Vj I-rn r) • 
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We claim that, for 1 < i j < m — r, g h] = g 3ti mod (gk,e, k > m — r), which implies 
that f re d verifies the statement. Let a it j denote the (i, j) — th entry of A(x). Let i < j and 
write 

m m 

9i,j .j T ^ cii/Uej and flj.i ®j,i T ^ ^ ®• 

l=m— r+1 l=m— r+1 

We deduce that g t) j — g 3 .i = J2T= m-r+i a i/Vt,j ~ a j/Vt,i since A is symmetric. Also, modulo 
the ideal (gk/, k > m — r), and for t > m — r + 1, one can explicit and a 3 j, by using 
polynomial relations gi ti = 0 and ggj = 0, as follows: 

m / m m 

9i,j 9j,i = ^ ^ ^ ^ a, i9 ty t ,mj ^ ^ cL ijt y t ,jy£,i 

£=m-r-\-l \ t=m—r +1 t=m—r +1 

m 

= ^2 a + ytjVtii) = 0 1110(1 (9k/, k > m — r). 

£,t=m-r -\-1 

The previous congruence concludes the proof. □ 

We prove below in Proposition [7] and in Corollary [TT] that, up to genericity assumptions, 
the ideal (/) = ( f re d ) is radical and that jj f re d matches the codimension of V r . In the next 
example, we explicitly write down the redundancies shown in Lemma [5] for a simple case. 



Example 6 We consider a 3x3 symmetric matrix of unknowns, and the kernel corre¬ 
sponding to the configuration {1,2} C {1,2,3}. Let 


( fll fl2 
/21 /22 
\ /31 /32 


X\ X 2 X 3 \ / 1 0 \ 

x 2 X A x 5 \ \ 0 1 

^3 x 5 x 6 J \ y 31 y 32 J 


We consider the classes of polynomials / 12,/21 modulo (f 3 i,f 3 2 ), deducing the following 
linear relation: f 12 - /21 = 2 / 32^3 - 2 / 31^5 = 2 / 31 ^ 2/32 - 2 / 32 ^ 2/31 = 0. 


Lagrange systems 

Let f(A , 1 ) be the polynomial system defining V r {A, t). We set c = m(m — r) + ) 

and e = ( m ~ r ), so that V r is defined by c = %f re d polynomial equations, and e = j \f — c 
is the number of redundancies eliminated by Lemma [5] We define, for M E GL n (C), the 
polynomial system £ = £(A o M, l), given by the coordinates of the map 

£ • ^ ^c+e _ y (£yi-\-m(rn—r)-\-c-\-e 

( x,y,z ) 1 —> (f(Ao M,t),z T Df(Ao M,l) - (ef,0)) , 

where e 3 E Q n is the first column of the identity matrix I„. We also define Z(Ao M, 1 ) = 
Z(£(A o M, t)). When V r (A o M, 1 ) is smooth and equidimensional, Z(A o M, 1 ) encodes 
the critical points of the restriction of Ili(a:, y) = x\ to V r (A o M, t). 
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Output representation 


The output of the algorithm is a rational parametrization (q 0 ,..., q n+ \) C Z[t\ such that 
the finite set defined in (jT]) contains at least one point on the spectrahedron SP. 

3.2 Real root finding for symmetric low rank loci 

We describe the main subroutine LowRankSym, which is a variant for symmetric pencils of 
the algorithms in [3TI :3Bj 3S]. Its output is a finite set meeting each connected component 
of T> r nl". It takes advantage of the particular properties of the incidence varieties over 
a symmetric low rank locus, as highlighted by Lemma [5j 


Properties 

We define the following properties for a given A 6 §^ +1 (Q): 

Property We say that A satisfies Pi if, for all i C {1,..., m}, with jj i — m — r, the 
incidence variety V r {A, i) is either empty or smooth and equidimensional. 

Property P 2 . We say that A satisfies P 2 if, for all r, T> r has the expected dimension 
n — ( m_ 2 +1 )- Property P 2 holds generically in §^ +1 (Q), as shown by Lemma HI 

We also define the following property for a polynomial system / C Q[x] and a Zariski 
open set ^cC": 

Property Q. Suppose that / C Q[x] generates a radical ideal and that it defines an 
algebraic set of codimension c, and let & C C* be a Zariski open set. We say that / 
satisfies Q in if the rank of Df is c in Z({f )) D &. 


Formal description of LowRankSym 

The formal description of our algorithm is given next. It consists of a main algo¬ 
rithm which checks the genericity hypotheses, and of a recursive sub-algorithm called 
LowRankSymRec. 
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LowRankSym 

Input: A G §^ +1 (Q), encoded by the m(m + l)(n + l)/2 
entries of A 0 , Ai,..., A n , and an integer 1 < r < m — 1; 
Output: Either the empty list [ ], if and only if T> r D R n = 
0, or an error message stating that the genericity assump¬ 
tions are not satisfied, or a rational parametrization q = 
(qo,..., q n+ i) C Z[f], such that for every connected com¬ 
ponent C C V r fl M n , with C fl X> r -i = 0, there exists 
t* G Z(q n+ 1 ) DM with (qi(t*)/q 0 (t*),..., q n (t*)/q 0 (t*)) G C. 

Procedure LowRankSym (A, r) 

1 . if n < ( m " r+1 ) then 

• if dim X> r = — 1 then return [ ], else return(“the 
input is not generic”); 

2. for l C {1,..., m} with %i — m- r do 

• if lsReg((74, l)) = false then return(“the input is 
not generic”); 

3. return(LowRankSymRec(A,r)). 

Procedure LowRankSymRec(74, r) 

1. choose M G GL n (Q); 

2. q ■<—[]; for i C {1,..., m } with )\i = m — r do 

• q b G- lmage(Project(RatPar(£(74 o M, t))), M _1 ); 

• q <- Union(g,gJ; 

3. choose t G Q; A i — (^4o ~L tA\, A 2 , • • •, ^4n)j 

4. q' <— Lift(LowRankSymRec(A, r), t); 

5. return(Union(g, q')). 


The routines appearing in the previous algorithm are described next: 

• IsReg. Input: A G S^ +1 (Q),t C {1 Output: true if V r (A,i) is empty or 

smooth and equidimensional of codimension m(m — r ) + ( m_ 2 r+1 ), false otherwise. 

• Project. Input: A rational parametrization of £(Ao M,t) C Q[x, y,z]; Output: an 
error message if the projection of Z(Ao M, 1 ) fl {(r, y , z ) : rank A(Mx ) = r} on the 
x— space is not finite; otherwise a rational parametrization of this projection. 
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• RatPar. Input: £(AoM,l) C Q[x, y,z\] Output: a rational parametrization of 
£(AoM,l). 

• Image. Input: a rational parametrization of a set Z C Q[xi,... , apv] and a matrix 
M G GLtv(Q); Output: a rational parametrization of M _1 Z = {a; G C w : Mx G 
Z}. 

• Union. Input: two rational parametrizations encoding Z X ,Z 2 C Q[xi,..., opy]; Out¬ 
put: a rational parametrization of Z x U Z 2 . 

• Lift. Input: a rational parametrization of a set Z C Q[xi,..., xn], and t G C; 
Output: a rational parametrization of {(£, x) : a; G Z}. 

3.3 Main algorithm: description 

The input of SolveLMI is a linear matrix A G 8^ +1 (Q), and the algorithm makes use of 
algorithm LowRankSym described previously, as a subroutine. The formal description is 
the following. 


SolveLMI 

Input: A G S" +1 (Q), encoded by the m(m + l)(n + l)/2 
entries of A 0 , A \,..., A n ; 

Output: Either the empty list [ ], if and only if {a; G 
M n : A(x) y 0} is empty; or an error message stating that the 
genericity assumptions are not satisfied, or, otherwise, either 
a vector x* = (x*,..., x* ) such that A(x*) = 0, or a rational 
parametrization q = ( q 0 ,..., q n+ \) C Z[f], such that there ex¬ 
ists t* G Z(q n+1 )nRwithA(q 1 (t*)/qo(t*),..,,q n (t*)/qo(t*)) h 
0. 

Procedure SolveLMI(d) 

1. x* <— SolveLinear(d); if x* yf [ ] then return(x*); 

2 . for r from 1 to m — 1 do: 

• q <— LowRankSym(^4, r); 

• if q — “the input is not generic” then return((/); 

• if q ^ [] then b <— CheckLMI(d, q); 

• if b — true then return(g); 

3. return([ ], “the spectrahedron is empty”). 


The different subroutines of SolveLMI are described next: 
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• SolveLinear. Input: A £ S^ +1 (Q); Output the empty list if A(x) = 0 has no solution, 
otherwise x* such that A(x*) = 0; 

• CheckLMI. Input : A £ S^ +1 (Q) and a rational parametrization q C Z[t]; Output 
true if A(qi(t*)/q 0 (t*),.. ., q n (t*)/qo{t*)) y 0 is satisfied for some t* £ Z(q n+1 ) flK, 
false otherwise. 

3.4 Main algorithm: correctness 

We prove in Theorem [TUI that SolveLMI returns a correct output if genericity properties 
on input data and on random parameters chosen during its execution are satisfied; the 
proof relics on some preliminary results that are described before. The proofs of these 
intermediate results share some techniques with [3TJ135., 36] and are given in Section |4j 


Intermediate results 

The first result is a regularity theorem for the incidence varieties V r {A,i). We focus on 
property Pi for the input matrix A (cf. Section 13721) . 

Proposition 7 Let m, n, r £ N, with 0 < r < m — 1. 

1. There exists a non-empty Zariski-open set sZ C S” +1 (C) such that if A £ srf D 
>§m +1 (Q)? then A satisfies Pi; 

2. if A satisfies Pi, there exists a non-empty Zariski open set AT C C such that if 
t £ fl Q, then (A 0 + tAi, A 2 ,..., A n ) satisfies Pi. 

Next, we prove that the projection of Z(A oM, t) D{(a;, y, z ) : rank A(M x) = r} over the 
x— space is finite and that this set meets the critical points of the restriction of the map 
III: (x, y) —> X] to V r (A, l). 

Proposition 8 Let A £ S"j +1 (Q) satisfy Pi. Then there exists a non-empty Zariski open 
set ^#i C GL n (C) such that, if M £ ^jflM n ,n(Q), f or all tc{l,..., m}, with jji = m—r, 
the following holds: 

1. The system £(A o M, t) satisfies Q in {(a;, y, z) : rank A(Mx) = r}; 

2. the projection of Z(AoM, t)n{(x, y, z) : rank4(Ma;) = r } on the x—space is empty 
or finite; 

3. the projection of Z(A o M, i) n {(x,y,z) : rank A(Mx) = r} on (x,y) contains 
the set of critical points of the restriction of IIi: (x,y) —> X\ to V r {A o M, i) fl 
{(x,y) : rank A(Mx) — r}. 
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Finally, we show, after a generic linear change of variables x, closure properties of the 
projection maps 7p(x) = (xi,... ,xf) restricted to V r . Also, in order to compute sample 
points on the connected components of T> r fl M n not meeting T> r _ i, the next proposition 
shows that to do that it is sufficient to compute critical points on the incidence variety 

V r . 

Proposition 9 Let A G §^ +1 (Q) satisfy Pi, and let d = dimD r . There exists a non¬ 
empty Zariski open set ^#2 C GL n (C) such that if M G .yf 2 DM ni „(Q), for any connected 
component C C V r fl M n , the following holds: 

1. for i = 1 ,d, 7 Tj(M _1 C) is closed; further, for t G K. lying on the boundary of 
7Ti(M^ 1 C), then 7r 1 ” 1 (f) fl M~ l C is finite; 

2. lett lie on the boundary of : for x G such that rank A(M x) = 

r, there exists m}, with jj t — m — r, and (x, y) G V r (A o M, l), such that 

III (x,y) = t. 


Theorem of correctness 

Let A G §m"|"m(Q) be the input of SolveLMI. We say that hypothesis H holds if: 

• A and all parameters generated by SolveLMI belong to the Zariski open sets defined 
in Proposition [?1 [8] and [9l for all recursive steps of LowRankSym; 

• A satisfies Property P 2 . 

We can now state the correctness theorem for SolveLMI. 

Theorem 10 (Correctness of SolveLMI) Suppose that H holds. Let TZ = {1 G 1" : A (a;) 
0} be the spectrahedron associated to A. Then two alternatives hold: 

1. TZ = 0; hence the output of SolveLMI with input A is the empty list; 

2. TZ 7 ^ 0; hence the output of SolveLMI with input A is either a vector x* such that 
A{x*) =0, if it exists; or a rational parametrization q = (q 0 ,, q n+ i) C Z [t\ such 
that there exists t* G Z(q n+ i ) flM with: 

• A(q 1 (t*)/q 0 (t *),..., q n (t*)/qo(t*)) h 0 and 

• rank A(qi(t*)/q 0 (t*),... ,q n (t*)/q Q (t*))—r(A) (cf. Notation \Tfj. 

Proof : Suppose Afx) = 0 has a solution. Hence, at Step [Tj of SolveLMI, SolveLinear with 
input A returns a vector x* such that A{x*) = 0. We deduce that x* G 5? 7 ^ 0 and that 
the rank of A attains its minimum on 5? at x*. We deduce that, if A{x) = 0 has at least 
one solution, the algorithm returns a correct output. 

Suppose now that either 5? is empty, or r(A) > 1. We claim that LowRankSym is correct, 
in the following sense: with input (A, r), with A satisfying Pi, the output is a rational 
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parametrization whose solutions meet each connected component C of D r nR n such that 

C H £> r _i = 0. 

We assume for the moment this claim, and consider two possible alternatives: 

1. ^ = 0. Consequently, CheckLMI outputs false at each iteration of Step El in 
SolveLMI. Thus the output of SolveLMI is the empty list, and correctness follows. 

2. 5? 7^ 0. Denote by C a connected component of V r (A) Cl R n such that C D 5? ^ 0. 
By Theorem [21 we deduce that C C 5?, and that C D U r (A)-i = 0- Let q be the 
output of LowRankSym at Step [21 of SolveLMI. By our claim, q defines a finite set 
whose solutions meet C, hence 5?. Consequently, CheckLMI returns true at Step [21 
and hence the algorithm stops returning the correct output q. 

We end the proof by showing that LowRankSym is correct. This is straightforwardly 
implied by the correctness of the recursive subroutine LowRankSymRec, which is proved 
below by using induction on the number of variables n. 

For n < ( m_ 2 +1 )’ s i n ce H holds, then A satisfies P r . Hence T> r is empty, and LowRankSym 
returns the correct answer [ ] (the empty list). 

Let n > ( m_ 2 r+1 )’ an( l (A,r) be the input. The induction hypothesis implies that for 
any A G §^(Q) satisfying Pi, then LowRankSymRec with input (A,r) returns a rational 
parametrization of a finite set meeting each connected component C C XVnR n_1 such that 
C fl V r -i = 0, with V r = {x G IT 1 " 1 : rankA(x) < r}. Let C C X\.nM n be a connected 
component with C fl T> r _i = 0, and let M be the matrix chosen at Step [U Hence, since H 
holds, by Proposition [9] the set 7 Ti(M _1 C) is closed in R. There are two possible scenarios. 

First case. Suppose that 7r 1 (M _1 C) = R, let t G Q be the rational number chosen at 
Step [3] of LowRankSymRec, and let A = (A 0 + tAi, A 2 ,..., A n ) G S”(Q). We deduce 
that 7rp 1 (t) fl M _1 C ^ 0 is the union of some connected components of the algebraic 
set T> r = {x G M n_1 : rankH(x) < r} not meeting D r -i- Also, since A satisfies Pi, 
so does A o M; by Proposition [71 then A satisfies Pi. By the induction assumption, 
LowRankSymRec with input ( A , r) returns at least one point in each connected component 
C C not meeting Z>, r _i, hence one point in C by applying the subroutine Lift at 

Step 01 Correctness follows. 

Second case. Otherwise, 7 Ti(M _1 C) ^ R and, since it is a closed set, its boundary is 
non-empty. Let t belong to the boundary of 7Ti (M _1 C), and suppose w.l.o.g. that 
7 Ti(M _1 C) C [t, +oo). Hence t is the minimum of the restriction of the map 7Ti to 
M~ 1 C. By Proposition El the set 7rp 1 (t) fl M~ 1 C is non-empty and finite, and for all 
x G 7r 1 " 1 (t) fl M _1 C, rank A(Mx) = r (indeed, for x G M _1 C, then M x G C and hence 
M x ^ V r _ i fl R"). Fix x G 7rp 1 (t) fl M _1 C. By Proposition El there exists i and 
y G r ) suc h that (x,y) G V r (AoM, i). Also, by Proposition [3 the set V r (AoM, l) is 

smooth and equidimensional. One deduces that (x, y ) is a critical point of the restriction 
of n x : (x, y) —» X\ to V r (doM, t) and that there exists 0 such that (x, y , z) G Z(AoM, i). 
Hence, at Step [21 the routine LowRankSymRec outputs a rational parametrization q L , 
among whose solutions the vector x lies. □ 
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4 Proof of intermediate results 


4.1 Proof of Proposition [7] 

We prove Assertion 1 and 2 separately. 

Proof of Assertion 1: Suppose w.l.o.g. that M = I n . For t C {l,...,m}, with 
j \t = m — r, let f re d C Q[x, y] be the system defined in Lemma [5] We prove that there 
exists a non-empty Zariski open set sY b C §^ +1 (C) such that, if A G fl S^ +1 (Q), 
f re d generates a radical ideal and Z(f red ) is empty or equidimensional, of codimension 
jj f r ed = rn(m — r) + ■ We deduce that, for A G A satisfies P l5 and we conclude 

by dehning sA = A b sY b (non-empty and Zariski open). 

Suppose w.l.o.g. that i — {1,..., m — r}. We consider the map 

ip : £n+m(m-r) x gri+lVQ _ ^ ^yn(m-r )+( m ~ 2 7 ' +1 ) 

(X, y, A) I-> f r e d 

and, for a fixed A G §^ + 1 (C), ^s section map <^ 4 : —» C m ^ m_r ^ + ( 2+ ) defined 

by p>A(x,y) = <p(x,y,A). Remark that, for any A, Zfpjf) = V r (A, t). 

Suppose </J - 1 (0) = 0: this implies that, for all A G Z(fred) — V r (A, t) = 0, that 

is A satishes Pi for A E sY b = §^ + 1 (C). 

If <£> - 1 (0) 7 ^ 0, we prove below that 0 is a regular value of ip. We conclude that by Thom’s 
Weak Transversality Theorem [SHI Section 4.2] there exists a non-empty and Zariski open 
set A b C §” + 1 (C) such that if A G srf L fl §™ +1 (Q), 0 is a regular value of ip a- Hence, 
by applying the Jacobian criterion (cf. [16] Theorem 16.19]) to the polynomial system 
/red, we deduce that for A G sY b D §” +1 (Q), V r (A, 1 ) is smooth and equidimensional of 
codimension jj f red . 

Let Dip be the Jacobian matrix of ip: it contains the derivatives of polynomials in f red 
with respect to variables x,y,A. We denote by a^j the variable encoding the th 

entry of the matrix Ag, £ = 0,..., n. We isolate the columns of Dip corresponding to: 

• the derivatives with respect to variables {do,i,j '■ i < m — r or j < m — r}; 

• the derivatives with respect to variables y t ,j such that i G 1 . 

Let (x, y, A) G <^ - 1 (0), and consider the evaluation of Dip at (x, y, A). The above columns 
contain the following non-singular blocks: 

• the derivatives w.r.t. {ao,i,j : i < m — r or j < m — r} of the entries of A(x)Y(y) 
after the substitution Y[ G- I m _ r , that is I( m - r )(m+r+i)/ 2 ; 

• the derivatives w.r.t. {y i} j : i G <■} of polynomials in Y[ — I m _ r , that is I( m _ r ) 2 . 

Hence, the above columns define a maximal non-singular sub-matrix of Dip at (x,y,A), 
of size m(m — r) + ( m b) +1 ) = jj f re d (indeed, remark that the entries of Y L — I m _ r do not 
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depend on variables a 0 ,i,j). Since (x,y,A) E p 1 (0) is arbitrary, we deduce that 0 is a 
regular value of p , and we conclude. □ 

Proof of Assertion 2: Fix i C {1 with j }i — m — r. Since A satisfies Pi, 

V. r (A, C) is either empty or smooth and equidimensional of codimension m(m— r) + ( m_ J +1 ). 
Suppose first that V r = 0. Hence for all t 6 C, V r fl {aq — t — 0} = 0, and we conclude 
by defining A? — C. Otherwise, consider the restriction of the projection map Hi : 
(x,y) —> xi to V r (A, l). By Sard’s Lemma [601 Section 4.2], the set of critical values of 
the restriction of 7T! to V r (A, l) is included in a finite subset PL C C. We deduce that, for 
t G A? — C \ PL, then ( A 0 + tA 1 , A 2 ,..., A n ) satisfies Pi. □ 

Corollary 11 Let sA C §^ +1 (Q) be as in Proposition [?] and let A G , 2 /. Then for every 
tc{l,..., m} with jjt = m — r, the ideal (f re d) = (/) zs radical, and V r (A, <.) is a complete 
intersection of codimension $f re d- 

Proof : We recall from the proof of Assertion 1 of Theorem [7] that, for A G sA , the 
rank of the Jacobian matrix of f re d is jj f re d = m{m — r) + ( m_ ^ +1 ) a t every point of 
V r (A, t). By the Jacobian criterion [[161 Theorem 16.19], the ideal (f re d) is radical and 
the algebraic set Z(f re d) = V r (A, l) is smooth and equidimensional of codimension jj f re d- 
Hence I(V r (A, t)) = ( f red ) can be generated by a number of polynomials equal to the 
codimension of V r (A, t), and we conclude. □ 


4.2 Proof of Proposition |8] 

Local equations of V r (A, t) 


Suppose A is a (not necessarily symmetric) linear matrix. Let us give a local description 
of the algebraic sets V r and V r (c/. also [34] Section 5]). Consider first the locally 
closed set T> r \ T > r -1 = {1 6 C n : rankA(x) = r}. This is given by the union of sets 
T> r fl G C n : det Nfx) 7 ^ 0} where N runs over all r x r sub-matrices of A(x). Fix 
i C {1,..., m} with = m — r. Let N be the upper left r x r sub-matrix of A(x), and 
consider the corresponding block division of A: 


A 


N Q \ 
P T R ) 


(3) 


with P,Q G M rm _ r (Q) and R G M m _ r)m _ r (Q). Let Q[x, |/]det v be the local ring obtained 
by localizing Q[x,|/] at (det A). Let yb) (resp. be the matrix obtained by isolating 

the first r (resp. the last m — r) rows of Y(y). Hence, the local equations of V r in 
{( x,y ) : det N(x) 7 ^ 0} are given by: 

Y w + iV" 1 Qy (2) = 0, S(A)y ( 2) = 0, y l - I m _ r = 0, (4) 


where E (N) = R — P T N~ l Q is the Schur complement of N in A. This follows from the 
following straightforward equivalence holding in the local ring Q[x, y] det tv (c/. also [34] 
Lemma 13]): 


A(x)Y(y)= 0 iff 


f l r 0 \ f N- 1 0 

V ~P T L-r ) V 0 l m- 


N Q 
P T R 


Y(y) = 0 . 
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Intermediate lemma 


Let w E C" be a non-zero vector, and consider the projection map induced by w. 
n w : (xi ,..., x n , y ) H- Wl x x H-b w n x n . 

For A E sY (given by Proposition [7]), for all l as above, the critical points of the restriction 
of n u , to V r (A, i ) are encoded by the polynomial system (/, g, h ) where 

/ = /(A, 1 ), ( 9 ,ft) = ^( on w ) =zT (y D o)’ (5) 

and z — (zi,, z c+e , 1) is a vector of Lagrange multipliers. Indeed, equations induced by 
( g , h ) imply that the vector w is normal to the tangent space of V r at (x, y ). 

We prove an intermediate lemma towards Proposition [HJ 

Lemma 12 Let A E §^ +1 (Q) satisfy Pi. Then there exists a non-empty Zariski open set 
W C C n such that, if w EW A Q n , for all l C {1,..., m}, j \i = m — r, the following 

holds: 

1. the system ( f,g,h ) in ([5]) satisfies Q in {( x,y,z ) : rankA(x) = r}; 

2. the projection of Z(f, g, h) D {(x, y, z) : rank Afx) = r} on the x—space is empty or 
finite; 

3. the projection of Z(f, g,h) A{(x,y,z) : rankA(;r) = r} on ( x,y ) contains the set of 
critical points of the restriction ofA w to V r D {(x,y) : rankA(:r) = r}. 

Proof of Assertion 1: The strategy relies on applying Thom Weak Transversality 
Theorem and the Jacobian criterion, as in the proof of Proposition [3 The similar passages 
will be only sketched. 

We claim (and prove below) that, given a r x r sub-matrix N of A, there exists a 
non-empty Zariski open set C C n such that, for w E Wm, (f,g,h) satisfies Q in 
{( x,y,z ) : det N ^ 0}. We deduce Assertion 1 by defining W = C\ n ^n, where N runs 
over all r x r sub-matrices of A(x). 

Let U L E C ( m ~ r ) xrn be the boolean matrix such that U L Y(y) = Y), and let U L = {U, V ' ! \ U ,, (2) ) 
be the subdivision with E £f m ~ r ) xr and Uj 2 ' 1 E The third equation in 

(HI) equals UfY(y ) — I m -r = 0. From ([4]) we deduce the equality 

i m _ r = u^y {1) + t/J 2 ) r (2) = (uj; 2) - u[ 1 ) n- 1 q)y { - 2) 

and hence that both W 2 ) an d [jj 2> _ \j [ 1J N~ 1 P are non-singular in Q[x, y]detN- We deduce 
that the above local equations of V r are equivalent to 

y (1) + iv _ 1 gy (2) = o, e (n) = o, y (2) - {u[ 2) - u^N^py 1 = o 

in Q[x, y]det n- We collect the above equations in a system /, of length c + e. Hence, the 
Jacobian matrix of / is 

( Dx [E(7V)]j j 0(m_r) 2 xm(m_r) \ 

Ir(m—r) * J • 

0 I(m—r) 2 / 
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Since A satisfies Pi, the rank of Df is constant and equal to c if evaluated at (x,y) £ 
Z(f) = V r (A, l) D {(x, y ) : det N 0}. Similarly to (jSJ), we define 


0, h) = z T 


Df 

w T 0 


with z = (zi ,..., z c+e , 1). The polynomial relations = 0 implies that Z( m _ r ) 2 +i = 0, 
for i — 1,... ,m(m — r ), and hence they can be eliminated, together with the variables 
Z( m _ r ) 2 +i , i = l,, m(m — r). Hence, one can consider the equivalent equations (/, g, h) 
where the last m(m — r) variables z do not appear in g. 

Let us define the map p: C n + c + e + m ( m -D xC"-> Qn+c+e+m(m-r), (p( x ,y, z,w) = ( f,g,h ), 
and for w £ C n , its section map <p w : (x, y, z) H y p{x, y, z, w). 

If <y5“ 1 (0) = 0, we define Wn = C n and conclude. Otherwise, let ( x,y,z,w ) £ <^ -1 (0). Re¬ 
mark that polynomials in / just depend on (x, y), hence their contribution in Dtp(x, y, z, w ) 
is the block Df , whose rank is c, since (x,y) £ V r . Hence, we deduce that the row- 
rank of Dip at (x,y,z,w) is at most n + c + m(m — r ). Further, by isolating the 
columns corresponding to the derivatives with respect to x,y, to uq,..., w n , and to 
Z(rn- r ) 2 +i, i = 1,..., m{m — r), one obtains a (n + c + e + m(m — r)) x (2 n + 2 m(m — r)) 
sub-matrix of Dip of rank n + c + m(m — r). □ 

Proof of Assertion 2: From Assertion 1 we deduce that the locally closed set £ = 
Z{f, g, h)fl{(;r, y, z) : rank A(x) = r} is empty or e—equidimensional. If it is empty, we are 
done. Suppose that it is e—equidimensional. Consider the projection map n x (x, y, z) = x, 
and its restriction to £. Let x* £ n x (£). Then rankA(a;*) = r and there exists a unique 
y £ C m(m " r) such that f(x*,y ) = 0. Hence the fiber 7T“ 1 (a:*) is isomorphic to the linear 
space {(zi,..., z c+e ) : ( Z \,..., z c+e )Df = {w T , 0)}. Since the rank of Df is c, one deduces 
that 7r' 1 (x*) is a linear space of dimension e, and by the Theorem on the Dimension of 
Fibers [SSI Sect. 6.3, Theorem 7] we deduce that 7r x (£) has dimension 0. □ 

Proof of Assertion 3: Since V r Fl {(x,y) : rankA(x) = r} is smooth and equidimen¬ 
sional, by [60 ; , Lemma 3.2.1], for w 0, the set crit (n w ,V r ) coincides with the set of 
points (x, y ) £ V r such that the matrix 


D(f,U w ) 


Df 

DU W 


has a rank < c. In particular there exists z = (z±, ..., z c+e , z c+e+ i) 0, such that 
z T D(f, n m ) = 0. One can exclude that z c+e+1 = 0, since this implies that Df has a 
non-zero vector in the left kernel, which contradicts the fact that A satisfies Pi. Hence 
without loss of generality we deduce that z c+e+ 1 = 1, and we conclude. □ 


Proof of the proposition 

We can finally deduce the proof of Proposition [HJ 

Proof of Proposition O Define as the set of matrices M £ GL„(C) such that the 
first row of ilW 1 is contained in W (defined in Lemma fT21) . The proof of all assertions 
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follows from Lemma IT2l since, for M £ , one gets 

/ Df(AoM,i) \ _ ( Df(A,i)oM M 0 \ 

V ef 0 • ■ ■ 0 ) \ w T 0 • • • 0 ) V o I m(m _ r) ) ’ 1 j 

where w T is the first row of M _1 . Indeed, for z = (z±, ..., z c+e ), we deduce from the 
previous relation that the set of solutions to the equations 


o' 

z T Df(A , 1 ) = (w T ,0) 

(7) 

is the image of the set of solutions of 

f(Ao M,l) — 0, 

z T Df(AoM,i) = (el,0) 

(8) 


by the linear map 



This last fact is straightforward since from ([6]) we deduce that system (J8J) is equivalent to 
f(A o M,l) = 0, z T ( Df(A , i) o M) = ( w T , 0). Hence the three assertions of Proposition [8] 
are straightforwardly deduced by those of Lemma [12j □ 

4.3 Proof of Proposition [9] 

For the proof of Assertion 1 of Proposition [91 we need to recall some notation from [Ml 
Sec. 5]. Let Z C C n be an algebraic set of dimension d. Its equidimensional component of 
dimensionp, for 0 < p < d, is Q p (Z). We define 5A{Z') = Q.q{Z) U- • •UOd_i(2)UsingO (i Z, 
where we recall that sing V denotes the singular locus of an algebraic set V, and 

d 

c 1o (7Tj, Z ) = Q 0 (Z) U ■ • • U Qi-i(Z) U crit ( 7 q, reg Q r Z), 

r=i 

Here reg V denotes V\sing V, 7q: x (aq ,,Xi) and crit (g, V) the set of critical points 
of the restriction of a map g to V. For M £ GL n (C) we recursively define 

• O d {M~ l Z) = M~ l Z ; 

• OiiM-'Z) = ^(Oi+^M-'Z)) U^( 7 r i+1 ,M- 1 Z) for i = 

0,..., d — 1. 

In [Ml Prop. 17] we proved that when M is generic in GL n (C) (that is, it lies out of 
a proper algebraic set) the algebraic sets Oi(M~ l Z) have dimension at most i and are 
in Noether position with respect to aq ,,Xi ( cf. [63], [16] for a background in Noether 
position). We used this fact in pH Prop. 18] to prove closure properties of the restriction 
of maps ni, i = 1,,d, to the connected components of Z fl M n . 
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Proof of Assertion 1: We denote by C GL n (C) the non-empty Zariski open set 
defined in [3U Prop. 17], for the algebraic set V r . Hence, for M G .#2 , we deduce by 
[341 Prop. 18] that for i = 1,..., d, and for any connected component C C T> r fl M n , the 
boundary of 7Tj(M _1 C) is contained in 7r ,((9j_i(M _1 D r ) C 7r i (M^ 1 C), and hence 

that 7 Tj(M _1 C) is closed. Moreover, let C C T> r fl M n be a connected component and let 
t G R be in the boundary of 7r 1 (M _1 C). Then [3U Lemma 19] implies that 7rf 1 (t) fl M _1 C 
is finite. □ 

Proof of Assertion 2: Let M G ^ 2 - Consider the open set 

6 = {(a;, 2 /) G C n+m( - m-r ) : rank A(M x) = r, rank Y(y) = m — r}. 

Let U x : (£ n + m { m - r ) n. c (a;,|/) = x. Then n x (^) is the locally closed set M _1 ('D r \ 

T> r _i) = {1 6 C" : rank A(M x ) = r}. We consider the restriction of polynomial equations 
A(M x)Y(y ) = 0 to 6. By definition of 6, we can split the locally closed set & fl 
Z(A(M x)Y(y )) into the union 

0nZ{A{Mx)Y(y)) = U (0 t nZ(A(Mx)Y(y))), 

i C {l,...,m} 

= m — r 

where ^ = {(x,y) : det YJ. 7 ^ 0 }. 

Let C be a connected component of T> r D M n . Let t lie in the frontier of 7Ti(M _1 C), and 
x G vr 7 1 (t) D M _1 C with rank A(M x) = r. Hence there exists t C {1,... ,m} such that 
x lies in the projection of V r (A o M, 1 ) on the x— space. Hence there exists y such that 
(x, y) G V r (A o M, l ) and such that y) = t. □ 


5 Complexity analysis 

Our next goal is to estimate the arithmetic complexity of algorithm SolveLMI, that is 
the number of arithmetic operations performed over Q. This will essentially rely on the 
complexities of state-of-the-art algorithms computing rational parametrizations. 

5.1 Output degree estimates 

We start by computing Multilinear Bezout bounds (c/. [60], Ch. 11]) on the output degree. 


Proposition 13 Let A G be the input of SolveLMI, and 0 < r < m — 1. Let 

p r = (m — r)(m + r + l)/2. If H holds, for all 1 C {1, ..., m}, with (Jt = m — r, the degree 
of the parametrization q L computed at Step [H of LowRankSymRec is at most 


9(m, 7i, r) 



n — 1 
— 1 — r{m 





5 
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with Gm : n, r = {k : max{0, n — p r } < k < iiiin{n — ( m 9 r+1 ), r(m — r)}}. Moreover, for all 
m,n,r, 9(m,n,r ) is bounded above by ( Pr ^ n ) 3 ■ 

Proof : We can simplify the system / = f(A, i) to a system of p r bilinear equations with 
respect to variables x = (aq,... ,x n ) and y = (y m - r+ yi,... ,ym,m-r )• Indeed, by Lemma 
El V r (H, l) is dehned by Y b — I m _ r = 0 and by m{m — r) — e — p r entries of A(x)Y(y), 
where e = ( m f r ) ■ Hence we can eliminate equations Y L — I m _ r = 0 and the corresponding 
variables y b j. Consequently, the Lagrange system can be also simplified, by admitting 
only p r Lagrange multipliers z. We can also eliminate the first Lagrange multiplier Z\ 
(since z =/=■ 0, one can assume z\ = 1) and impose a rank defect on the truncated Jacobian 
matrix obtained by Df by eliminating the first column (that containing the derivatives 
with respect to aq). 

The bound 9(m,n,r), by (601 Ch. 11], equals the coefficient of s p z r ~ l in ( s x + 

s y) Pr ( s y + s z) n-1 ( s ;r + s z )A m - r \ This can be easily obtained by writing down such an 
expansion and solving the associated linear system forcing the constraints on the expo¬ 
nents of the monomials. The result is exactly the claimed closed formula. The estimate 
9{m,n,r) < ( Pr + n ) is obtained by applying the standard formula: 



with a — n and b = p r . □ 

Corollary 14 With the hypotheses and notations of Proposition [I3J the sum of the de¬ 
grees of the rational parametrizations computed during SolveLMI is bounded above by 
J2 r <r(A) (7)^( m ’ n 5 r )- The degree of the rational parametrization whose solutions intersect 
5^ is in 



Proof : Remark that the number of subsets i C {1,..., m}, with §l = m — r is ( n ' l 'f r ) = 
(™), and that SolveLMI stops when r reaches r(A). Hence the first part follows by applying 
Proposition E31 Finally, remark that Po > Pi > ■ ■ ■ > p r > ■ ■ ■ for all m, and hence, by 
Proposition E31 the degree of the rational parametrization whose solutions intersect 5? is 

of the order of ( Pr + n ) 3 < ( P0 + n ) 3 = ■ □ 

In the column deg of Table |T] we report the degrees of the rational parametrization q L 
returned by LowRankSymRec at Step [2] compared with its bound 9{m,n,r) computed in 
Proposition E2J For this table, the input are randomly generated symmetric pencils with 
rational coefficients. When the algorithm does not compute critical points (that is, when 
the Lagrange system generates the empty set) we put deg = 0. 

We recall that the routine LowRankSymRec computes points in components of the real 
algebraic set T> r D M n not meeting the subset V r -1 FI M n , hence of the expected rank r. 
Moreover, we recall that LowRankSym calls recursively its subroutine LowRankSymRec, 
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eliminating at each call the first variable. Hence, the total number of critical points 
computed by LowRankSym for a given expected rank r is obtained by summing up the 
integer in column deg for every admissible value of n. We remark here that both the degree 
and the bound are constant and equal to 0 if n is large enough. Hence, the previous sum is 
constant for large values of n. Similar behaviors appear, for example, when computing the 
Euclidean Distance degree (EDdegree) of determinantal varieties, as in [T5j or [50] . In [501 
Tablet], the authors report on the EDdegree of determinantal hypersurfaces generated 
by linear matrices A(x) = A 0 + X\Ai + • • • + x n A n : for generic weights in the distance 
function, and when the codimension of the vector space generated by Hi,..., A n is small 
(for us, when n is big, since matrices Ai are randomly generated, hence independent for 
n < = clini8 m (Q)) the EDdegree is constant. Similar comparisons can be done 

with data in [50], Example 4] and [50] Corollary 3.5]. 


(m, r, n) 

deg 

9(m, n, r) 

(m, r, n ) 

deg 

6{m, 7i, r) 

(3,2,2) 

6 

9 

(4,3,9) 

0 

0 

(3,2,3) 

4 

16 

(5,2,5) 

0 

0 

(3,2,4) 

0 

15 

(5,2,6) 

35 

924 

(3,2,5) 

0 

6 

(5,2,7) 

140 

10296 

(3,2,6) 

0 

0 

(5,3,3) 

20 

84 

(4,2,3) 

10 

35 

(5,3,4) 

90 

882 

(4,2,4) 

30 

245 

(5,4,2) 

20 

30 

(4,2,5) 

42 

896 

(5,4,3) 

40 

120 

(4,2,6) 

30 

2100 

(5,4,4) 

40 

325 

(4,2,7) 

10 

3340 

(5,4,5) 

16 

606 

(4,2,8) 

0 

3619 

(6,3,3) 

0 

0 

(4,2,9) 

0 

2576 

(6,3,4) 

0 

0 

(4,2,12) 

0 

0 

(6,3,5) 

0 

0 

(4,3,3) 

16 

52 

(6,3,6) 

112 

5005 

(4,3,4) 

8 

95 

(6,4,2) 

0 

0 

(4,3,7) 

0 

20 

(6,4,3) 

35 

165 

(4,3,8) 

0 

0 

(6,5,3) 

80 

230 


Table 1: Degrees and bounds for rational parametrizations 


The values in column deg of Table |T] must also be compared with the associated algebraic 
degree of semidefinite programming. Given integers k,m,r with r < m — 1, Nie, Ranes- 
tad, Sturmfels and von Bothmer computed in [481 [24 ] formulas for the algebraic degree 
5{k, m , r) of a generic semidefinite program associated to mx m k— variate linear matrices, 
with expected rank r. Since the values in column deg match exactly the corresponding 
values in [48] Table 2], we conclude this section with the following expected result, which 
is a work in progress. 

Conjecture 15 Let A e §^ +1 (Q) be the 'input o/SolveLMI, and suppose that y = {i6 
M n : A(x) h 0} is not empty. Let 5(k,m,r ) be the algebraic degree of a generic semidefi¬ 
nite program with parameters k,m,r as in If H holds, then the sum of the degrees 
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of the rational parametrizations computed during SolveLMI is 

r(A) , v min (n,p r -\-r(m—r)) 

Y ( 7 ) Y S(k,m,r), 

r= 1 k=p r -r(m-r) 

where p r = (m — r)(m + r + l)/2. 

5.2 The complexity of SolveLMI 

Complexity of some subroutines 

We first provide complexity estimates for subroutines SolveLinear, CheckLMI, Project. Lift, 
Image and Union. 

Complexity of SolveLinear. This subroutine computes, if it exists, a solution of the A(x) = 
0. This can be essentially performed by Gaussian elimination. The complexity of solving 
( m ^ 1 ) linear equations in n variables is linear in ( m + 1 ) and cubic in n. 

Complexity of CheckLMI. This subroutine can be performed as follows. Let q C Z[t\ 
be the rational parametrization in the input of CheckLMI. The spectrahedron y = {x G 
M n : A(x) 0} is the semi-algebraic set defined, e.g., by sign conditions on the coefficients 
of the characteristic polynomial 

p(s,x) = det(H(x) + sl m ) = f m (x) + f m -i(x)s H-f /i(x)s m_1 + s m . 

That is, y = {x G R" : /j(x) > 0,Vi = 1 We make the substitution x* «— 

qi(t) / qo(t) in A{x) and compute the coefficients of p(s,x(t)), that are rational functions 
of the variable t. Hence CheckLMI boils down to deciding on the sign of m univariate 
rational functions (that is, of 2m univariate polynomials) over the finite set defined by 
q n+ i(t) = 0. We deduce that the complexity of CheckLMI is polynomial in m and on the 
degree of q n+1 (that is, on the degree of q ) see [6| Ch. 13]. 

Complexity of Project, Lift, Image and Union Estimates for the arithmetic complexities 
of these routines are given in [60], Ch. 10]. In particular, if 6 = 6(m,n,r) is the bound 
computed in Proposition [13l and h — n + r(m — r) + p r , then: 

• By Lemma 10.1.5], Project is performed within 0~{n 2 6 2 ) operations; 

• By [60, Lemma 10.1.6], Lift is performed within O~(h0 2 ) operations; 

• By [60^ Lemma 10.1.1], Image is performed within 0~(n 2 9 + h 3 ) operations; 

• By 00. Lemma 10.1.3], Union is performed within 0~{n6 2 ) operations. 


Complexity of LowRankSym and SolveLMI 

The complexity of LowRankSym boils essentially down to the complexity of LowRankSym- 
Rec, that is the complexity of RatPar. This is performed with the symbolic homotopy 
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algorithm [39]: we bound its complexity in this section. We just remark that computing 
the dimension of T> r at Step [Tj of LowRankSym and performing the control routine IsReg 
can be done by combining the Jacobian criterion and Grobner bases computations. Onr 
complexity analysis omits this step. 

We recall that the simplified Lagrange system defined in the proof of Proposition [T31 
contains: p r = (m — r)(in + r + l)/2 polynomials of multidegree bounded by (1,1,0); 
n— 1 polynomials of multidegree bounded by (0,1,1); r(m — r ) polynomials of multidegree 
bounded by (1, 0,1). Let us denote by £ this system. We denote by 

A xy = {1 pXi.yj.XiVj : i = 1,... ,n,j = 1,..., r(m - r)} 

A yz = {l,yj,z k ,yjZ k : j = 1,..., r(m - r), k = 2 
A xz = {1 ,Xi,z k ,XiZ k : i = 1,.. .,n, k = 2,... ,p r } 

the supports of the aforementioned three groups of polynomials. Let £ C Q[x, y, z] be a 
polynomial system such that: 

• the length of £ equals that of t\ 

• for i — 1,..., n — 1 + m 2 — r 2 , the support of G equals that of £p 

• the solutions of £ are known. 

Remark that £ can be easily built by considering suitable products of linear forms. We 
build the homotopy 

t£+ (1 - t)£ C Q[x,y,z,r], (9) 

where r is a new variable. The system (J9J) defines a curve. From [391 Proposition 6.1], if the 
solutions of £ are known, one can compute a rational parametrization of Z(t£ + (1 — t)£) 
within 0{{n 2 N\ogQ + n U}+1 )ee') arithmetic operations over Q, where: n is the number 
of variables in £\ N = p r jjA xy + (n — l)#A yz + r(m — r)$A xz ] Q = max{||g|| : q e A xy U 
A yz U A xz }; e is the number of isolated solutions of i\ e' is the degree of Z(t£+ (1 — t)£)\ 
uj is the exponent of matrix multiplication. 

The technical lemma below gives a bound on the degree of Z{t£ + (1 — 1)£). 

Lemma 16 Let Q m ,n,r and 0(m,n,r ) be the set and the bound defined in Proposition 
m and suppose that Gm, n ,r 0- Let e! be the degree of Z(t£+ (1 — t)£). Then e' G 
O (fin + p r + r(m — r)) min{n, p r ) 8(m, n, r)). 

Proof : Similarly to Proposition [131 we exploit the multilinear structure of t£ + (1 — t)£, 
to compute a bound on e'. The system is bilinear with respect to x,y,z,r. We recall 
jjx = nfiy = r(m — r)fiz = p r — 1, jjr = 1, with p r = (m — r)(m + r + l)/2. By ^9, 
Ch. 11], e' is bounded by the sum of the coefficients of q — (s^ + s y + s T ) Pr (s y + s z + 
Sr) n_1 (sx + S z + s T ) r(m_T ' ) modulo I = (s£ +1 , s W' t_r ) +1 , s Pr^ s 2 ) c z[s x ,Sy,s z ,s T ). We see 
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that q = qi + s r (q 2 + 93 + 94 ) + g with s 2 that divides g and 


9i = (s x + s y r(s y + s z ) n -\s x + s z ) r(m " r) 

g 2 = p r (s x + s,,)*- 1 ^ + s z ) n ~\ Sx + Sz y (w *- r) 

Qs = ( U ~ l)(Sx + Sy) Pr (Sy + S^) n 2 (S X + S,) r(m ^ 

94 = r(m - r)( s * + s t/) Pr ( s y + Sz) n_1 (sx + s z y( m - r) -\ 

Hence q = qi + s r (g 2 + 93 + < 74 ) mod /, and the bound is given by the sum of the 
contributions of 91 , 92,93 and 94 (multiplying 92 , 93,94 by s T does not change the sum of 
the coefficients modulo I). The contribution of 91 is the sum of the coefficients of its 
class modulo I' = (s” +1 , s r y m r ^ +1 ,s Pr ). This has been computed in Proposition fl3l and 
coincides with 9(m, n, r ). We compute the contribution of q 2 . Let q 2 = p r q 2 • It is sufficient 
to compute the sum of the coefficients of q 2 modulo I' (defined above), multiplied by p r . 
Since deg 92 = n — 2 + p r + r(m — r), and since the maximal powers admissible modulo 
/' are and s Pr ~ l , three configurations are possible. 

(A) The coefficient of s ™ _1 s T y ( ' m ~ r ' > sp ' -1 in q 2 , that is 

s = vf Pr_1 V n ~ l \ / r(m — r)\ 

A \n — 1 — k) \k — 1 + p r — r(m — r)J \ k )' 

(B) the coefficient of s™s r y m r) 1 s Pr ~ l in q 2 , that is 

^ > = spfPr~ 1\ / n- 1 \fr(m-r )\ 

^ \n — k) \k — 1 + p r — r(m — r)J \ k J ’ 

(C) the coefficient of s Pr ~ 2 in q 2 , that is 

% = ^fPr~ 1\ / n- 1 \/r(m-r)\ 

y^\n-kJ\k-2+p r - r(m -r)J\ k J 

Hence we need to bound the expression p r (J^A + E# + Sc). One can easily check that 
Sa < 9(m,n,r) and E# < 6(m,n,r), while the same inequality is false for Ec- However, 
we claim that Ec < (1 + min {n,p r }) 9(m,n : r) and hence that the contribution of q 2 is 
Pri^A + E b + Ec) G O ( p r min{n,p r } 6 ( 771 , n, r)). We prove now this claim. Let 

Xi = max{0, n — p r } X 2 — minjn — p r + r(m — r),r(m — r)} 

ol\ = max{0, n — p r + 1} a 2 = min{n — p r + r(m — r ) + 1, r(m — r)} 

so that the sum in 9(m,n, r ) runs over xi < k < y 2 and that in Ec over aq < k < a 2 . 
Remark that Xi A aq and y 2 < a 2 . Denote by <p(k) the k— th term in the sum dehning 
Ec, and by y(fc) the k— th term in the sum dehning 9(m,n,r). Then for all indices k, 
admissible both for 9(m , n, r ) and Ec, that is for aq < k < y 2 , one gets, by basic properties 
of binomial coefficients, that tp(k) = ^/(k)x(k), with T(fc) = • When k 



runs over all admissible indices, T(fc) is non-decreasing monotone, and its maximum is 
tf(* 2 ) and is bounded by min{n,p r }. We deduce the claimed inequality Y>c < (1 + 
min{n,p r }) 9{m,n,r), since if y 2 < ay then y 2 = a 2 — 1 and p(a 2 ) is bounded above by 
9{m, n, r). 

Contributions of q 3 and q A . As for q 2 , we deduce that the contribution of q 3 is in 
O (n inin{n, p r } 9(m, n, r )) and that of q 4 is in O ( r{m — r ) min{n,p r } 9{in, n , r)). □ 

We use this degree estimate to conclude our complexity analysis of LowRankSym. 

Proposition 17 Let A e §^ +1 (Q) be the input of SolveLMI and 0 < r < m — 1. Let 
9[m,n,r ) be the bound defined in Proposition [731 Let p r = (m — r)(in + r + l)/2. Then 
RatPar returns a r.p. within 0~ (( m ) (n + p r + r(m — r )) 7 9(m,n,r) 2 ) arithmetic opera¬ 
tions over Q. 


Proof : Let l be the simplified Lagrange system as in the proof of Proposition [13j We 
consider the bound on the degree of the homotopy curve given by Lemma [16j We deduce 
the claimed complexity result by applying [391 Proposition 6.1], and by recalling that 
there are many subsets of { 1 ,..., m} of cardinality m — r. □ 

We straightforwardly deduce the following complexity estimate for SolveLMI. 


Theorem 18 (Complexity of SolveLMI) Let A 6 §^ + 1 (Q) be the input symmetric pen¬ 
cil and suppose that H holds. Then the number of arithmetic operations performed by 

SolveLMI are in 


O' 





{n + (m — r){m + 3r)) 


Proof : From Proposition [T71 we deduce that LowRankSymRec runs essentially within 
0 '((”) (n+p r +r(m — r )) 7 9{m,n,r ) 2 ) arithmetic operations. The inequality 9(m,n,r ) < 
(n+p r ) 3 - g p rovec | j n Proposition [13j Moreover, there are at most n recursive calls of 
LowRankSymRec in LowRankSym, and SolveLMI stops at most when r reaches m — 1. Fi¬ 
nally, the cost of subroutines SolveLinear, CheckLML Project, Lift, Image and Union is negli¬ 
gible. We deduce that the complexity of SolveLMI is in O' ^n X^ r <m-i (7) ( n + Pr + r ( m ~ r 
Since p r < Po = m 2 f m and p r + r(m — r) < (m — r){m + 3r), we conclude. □ 


6 Experiments 

SolveLMI is implemented in a maple function, and it is part of a library called spectra 
(Semidefinite Programming solved Exactly with Computational Tools of Real Algebra), 
released in September 2015. It collects efficient and exact algorithms solving a large class 
of problems in real algebraic geometry and semidefinite optimization. 

We present in this section the results of our experiments. These have been performed 
on a machine with the following characteristics: Intel(R) Xeon(R) CPU E7540@2.00GHz 
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with 256 Gb of RAM. We rely on the maple implementation of the software FGb mi. 
for fast computation of Grobner bases. To compute the rational parametrizations we use 
the implementation in MAPLE of the change-of-ordering algorithm FGLM [2TJj| and of its 
improved versions EDug. 

6.1 Generic symmetric pencils 

We implemented the function LowRankSym and tested the running time of the imple¬ 
mentation with input generic symmetric linear matrices. We recall that the algorithm 
SolveLMI amounts to iterating LowRankSym by increasing the expected rank r. By generic 
data we mean that a natural number N e N large enough is fixed, and numerators and 
denominators of the entries of A^t = 0are uniformly generated in the interval 
[-N, TV], We report in Table [2] the timings and the degrees of output rational parametriza¬ 
tions. 


(■ m , r, n) 

PPC 

LRS 

totaldeg 

deg 

(to, r, n) 

PPC 

LRS 

totaldeg 

deg 

(3,2,2) 

0.2 

8 

9 

6 

(4,3,9) 

OO 

28 

40 

0 

(3,2,3) 

0.3 

11 

13 

4 

(4,3,10) 

OO 

29 

40 

0 

(3,2,4) 

0.9 

13 

13 

0 

(4,3,11) 

oo 

30 

40 

0 

(3,2,5) 

5.1 

14 

13 

0 

(5,2,2) 

0.6 

0 

0 

0 

(3,2,6) 

15.5 

15 

13 

0 

(5,2,3) 

0.9 

0 

0 

0 

(3,2,7) 

31 

16 

13 

0 

(5,2,4) 
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1 

0 

0 

(3,2,8) 

109 

17 

13 

0 

(5,2,5) 

1.6 

1 

0 

0 

(3,2,9) 
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18 

13 

0 

(5,2,7) 

oo 

25856 
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(4,2,2) 

0.2 

0 

0 

0 

(5,3,2) 

0.4 

1 

0 

0 

(4,2,3) 

0.3 

2 

10 

10 

(5,3,3) 

0.5 

3 

20 

20 

(4,2,4) 

2.2 

9 

40 

30 

(5,3,4) 

OO 

1592 
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90 

(4,2,5) 

12.2 

29 

82 

42 

(5,3,5) 

oo 

16809 

317 

207 

(4,2,6) 

OO 

71 

112 

30 

(5,4,2) 

0.5 

7 

25 

20 

(4,2,7) 

oo 

103 

122 

10 

(5,4,3) 

10 

42 

65 

40 

(4,2,8) 

oo 

106 

122 

0 

(5,4,4) 

OO 

42 

105 

40 

(4,2,9) 

oo 

106 

122 

0 

(5,4,5) 

oo 

858 

121 

16 

(4,3,3) 

1 

10 

32 

16 

(6,3,3) 

4 

0 

0 

0 

(4,3,4) 
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21 

40 

8 

(6,3,4) 

140 

1 

0 

0 

(4,3,5) 

oo 

22 

40 

0 

(6,3,5) 

oo 

2 

0 

0 

(4,3,6) 

oo 

24 

40 

0 

(6,3,6) 

oo 

704 

112 

112 

(4,3,7) 

oo 

26 

40 

0 

(6,4,2) 

0.6 

1 

0 

0 

(4,3,8) 

oo 

27 

40 

0 

(6,5,3) 

OO 

591 

116 

80 


Table 2: Timings and degrees for dense symmetric linear matrices 


We recall that m is the size of the input matrix, n is the number of variables and r is 
the expected maximum rank (that is, the index of the algebraic set T> r ). We compare our 
timings (reported in LRS) with those of the function PointsPerComponents (column PPC) 
of the library RAGLIB developed by the third author [56]. The input of PointsPerCompo¬ 
nents are the (r + 1) x (r + 1) minors of the linear matrix, and the output is a rational 
parametrization of a finite set meeting each connected component of T> r n M n . We do 
not consider the time needed to compute all the minors of A{x) in PPC. The symbol oo 
means that we did not succeed in computing the parametrizations after 48 hours. Column 
deg contains the degree of the parametrization returned by LowRankSymRec at Step [2j 
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or 0 if the empty list is returned. Column totaldeg contains the sum of the values in deg 
for k varying between 1 and n. For example, for m — 4, r = 2, for n < 2 and n > 8 
the algorithm does not compute critical points, while it computes rational parametriza- 
tions of degree respectively 10, 30,42, 30,10 for n — 3,4, 5, 6 , 7; the number 82 in totaldeg 
for (m, n, r ) = (4, 2, 5) is obtained as the sum 10 + 30 + 42 of the integers in deg for 
77i = 4, r = 2 and n = 3,4,5. We remark that, as for Table [T] the value in deg for a 
given triple m, n, r coincides with the algebraic degree of semidehnite programming, that 
is 5(n,m,r ) as defined in [48]. 

Our algorithm allows to tackle examples that are out of reach for raglib and that, most 
of the time, the growth in terms of running time is controlled when m, r are fixed. This 
shows that our dedicated algorithm leads to practical remarkable improvements: indeed, 
for example, 4x4 linear matrices of expected rank 2 are treated in a few minutes, up to 
linear sections of dimension 9; we are also able to sample hypersurfaces in M 5 defined by 
the determinant of 5 x 5 symmetric linear matrices; finally, symmetric linear matrices of 
size up to 6 with many rank defects are shown to be tractable by our approach. We observe 
that most of the time is spent to compute a Grobner basis of the Lagrange systems, and 
for this we use new fast algorithms for the change of monomial orderings |19j . 

6.2 Scheiderer’s spectrahedron 

We consider the symmetric pencil 


/ 

1 

0 

Xl 

0 

—3/2 — X2 

X3 

\ 


0 

— 2x\ 

1/2 

X2 

1 

to 

1 

4^ 

-X5 



Xi 

1/2 

1 

X4 

0 

Xq 



0 

X2 

£4 

2*3 + 2 

x 5 

1/2 



—3/2 — x 2 

—2 — X4 

0 

x$ 

2^0 

1/2 


V 

X3 

-X5 

x 6 

1/2 

1/2 

1 

/ 


which is the Gram matrix of the trivariate polynomial 

f(ui, U2, U3) = u\ + U\u\ + U2 — 3ulll2U3 — Au\u\u3 + 2 ulul + U\u\ + U2u\ + u\. 

In other words, / = v T A(x)v for all x G M 6 , where v = (u‘l,U\U 2 ,U 2 ,UiU 3 ,U 2 U 3 ,uf) 
is the monomial basis of the vector space of homogeneous polynomials of degree 2 in 
Ui,u 2 ,u 3 . Since / is globally nonnegative, by Hilbert’s theorem [[38] it is a sum of at 
most three squares in M[iti, it 2 , w 3 ], namely there exist f 2 , h £ ®[wi,W 2 ,w 3 ] such that 
/ = f f + /| + /[|. Moreover, the spectrahedron y = {1 6 K 6 : i(i) ^ 0} parametrizes 
all the sum-of-squares decompositions of / (and it is a particular example of a Gram 
spectrahedron). Hence 5? must contain a matrix of rank at most 3. 

Scheiderer proved in [ 6 T] that / does not admit a sum-of-squares decomposition in the ring 
Q[iti, U 2 , U 3 ], that is, the summands / 1 , / 2 , / 3 cannot be chosen to have rational coefficients, 
answering a question of Sturmfels. We deduce that S? does not contain points with 
rational coordinates. In particular, it is not full-dimensional (its affine hull has dimension 
< 5) by straightforward density arguments. 

We first easily check that V 0 fl M 6 = V\ D M 6 = 0). Further, for r = 2, the algorithm 
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returns the following rational parametrization of D 2 fl R 6 : 

/ 3 + 16 1 8 - 24 1 2 8 + 6 1 + 8 1 2 16 + 6 1 - 16 1 2 -3 - 16* 3 + 16t \ 

\ —8 + 24t 2 ’ —8 + 24i 2 ’ -8 + 24t 2 ’ -8 + 24t 2 ’ -8 + 24i 2 ’ -8 + 24i 2 ) ' 

where t satishes 8 1 3 — 8t — 1 = 0. The set X> 2 is, indeed, of dimension 0, degree 3, and it 
contains only real points. Remark that the technical assumption P 2 is not satisfied here, 
since the expected dimension of T> 2 is —1. Conversely, the regularity assumptions on the 
incidence varieties are satisfied. By applying CheckLMI one gets that two of the three 
points lie on that is those with the following floating point approximation up to 9 
certified digits: 


0.127050844 \ 

1.000000000 
0.967716166 
1.967716166 ' 

0.127050844 
0.127050844 

These correspond to the two distinct decompositions of / as a sum of 2 squares. An ap¬ 
proximation of such representations can be computed by factorizing the matrix A(x(t*)) = 
V T V where t* is the corresponding root of 8 1 3 — 8t — l and V e M 2 ,6(M) is full rank. The 
corresponding decomposition is / = v T V T Vv = ||Ru|| 2 . At the third point of V- 2 HR 6 the 
matrix A(x) is indefinite, so it is not a valid Gram matrix. 

To conclude, SolveLMI allows to design a computer-aided proof of Scheiderer’s results. 
This example is interesting since the interior of 5? is empty and, typically, this can lead 
to numerical problems when using interior-point algorithms. 


/ -0.930402926 \ 
- 1.000000000 
0.731299211 
-0.268700788 
0.930402926 
\ -0.930402926 J 


( ~ 


and 
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7 Conclusion 

We have presented a probabilistic exact algorithm that computes an algebraic represen¬ 
tation of at least one feasible point of a LMI A{x) cz 0, or that detects emptiness of 
5^ = {x E R n : A{x) cz 0}. The algorithm works under assumptions which are proved 
to be generically satished. When these assumptions are not satisfied, the algorithm may 
return a wrong answer or raises an error (when the dimension of some Lagrange system 
is not 0). The main strategy is to reduce the input problem to a sequence of real root 
finding problems for the loci of rank defects of A(x)\ if 5^ is not empty, we have shown 
that computing sample points on determinantal varieties is sufficient to sample 5 ?, and 
that it can be done efficiently. Indeed, the arithmetic complexity is essentially quadratic 
on a multilinear Bezout bound on the output degree. 

This is, to our knowledge, the first exact computer algebra algorithm tailored to linear 
matrix inequalities. We conjecture that our algorithm is optimal since the degree of the 
output parametrization matches the algebraic degree of a generic semidefmite program, 
with expected rank equal to the minimal achievable rank on 5?. Since deciding the 
emptiness of 5? is a particular instance of computing the minimizer of a linear function 
over this set (namely, of a constant), our algorithm is able to compute minimal-rank 
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solutions of special semidefinite programs, which is, in general, a hard computational task. 
Indeed, numerical interior-point algorithms typically return approximations of feasible 
matrices with maximal rank among the solutions (those lying in the relative interior of 
the optimal face). Moreover, the example of Scheiderer’s spectrahedron shows that we can 
also tackle degenerate situations with no interior point which are typically numerically 
troublesome. 

To conclude, as highlighted by the discussions in Section [6] our viewpoint includes an 
effective aspect, by which it is essential to translate into practice the complexity results 
that have been obtained. This is the objective of our MAPLE library SPECTRA. It must 
be understood as a starting point towards a systematic exact computer algebra approach 
to semidefinite programming and related questions. 
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